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In a quenched mesoscopic fluid, modelling transport processes at high densities, we perform 
computer simulations of the single particle energy autocorrelation function C e (t), which is essentially 
a return probability. This is done to test the predictions for power law tails, obtained from mode 
coupling theory. We study both off and on-lattice systems in one- and two-dimensions. The predicted 
long time tail ~ t~ d ^ 2 is in excellent agreement with the results of computer simulations. We also 
account for finite size effects, such that smaller systems are fully covered by the present theory as 
well. 
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I. INTRODUCTION 

In classical fluids in thermal equilibrium the velocity autocorrelation function (VACF) decays 
algebraically ~ t~ d / 2 at large times, and so do all Green-Kubo integrands, as established by 
molecular dynamics 0- 12- Q ■ mode coupling 0, [j| , and kinetic theory Q . These functions are 
equilibrium time correlation function (J(i)J(0)) o of single- or iV-particle currents, whose time 
integrals give the transport coefficients, such as self-diffusion, viscosity or heat conductivity. 

The long time tails are not restricted to current-current correlations, but apply to large classes 
of time correlations, such as single site correlations, provided the dynamics obeys a conservation 
law. Similar single particle correlations with a long time tail (LTT) exist in fluids p|. These tails 
have a rather universal shape, in the sense that they are determined by the decay of the slow 
macroscopic modes of the system, and are independent of the microscopic details. 

Recently, we have discussed in paper I a simplified mesoscopic model for transport in fluids, 
in which the rapid short range fluctuations are averaged out. It is referred to as random DPD solid 
(Dissipative Particle Dynamics), and defined by a fluctuating heat equation (Langevin equation), 
where the random force satisfies the fluctuation-dissipation theorem. It is a special case of general 
DPD-fluids 0, in which point particles are characterized by positions r,*, velocities v$ and 
possibly by an internal energy q [i = 1, 2, • • • , N). 

The simplification has been obtained by quenching the translational degrees of freedom of 
the DPD particles. The only dynamic degrees of freedom left are the internal energies €i(i — 
1, 2, • • • , N) of the particles, where the total energy £\ e^(t) = E is conserved. The only remaining 
transport process is heat diffusion between fixed particles, where energy hops instantaneously be- 
tween interacting pairs {ij} within the interaction range (r^ < r c ). This transport mechanism is 
called collisional transfer. The direction of the energy flow is determined by the "local temperature 
gradient" [cj — ej), and heat conduction is only possible for densities above a percolation threshold 
Pp. In the underlying percolating structure two particles are connected (by a bond) if the distance 
between them satisfies, < r c . So we are dealing with bond percolation diffusion on the random 
proximity network, where the transition rate or conductivity across a bond is constant. 

In this paper we study the energy auto-correlation function of a single particle C e (t). It is the 
analog of the probability of return to the origin, P(R = 0,t), where P(R,t) satisfies a diffusion 
equation. 

The random DPD model is complementary to the Lorentz gas with overlapping scatterers in 
several respects. The latter has only kinetic transport of particles and the former has only the 
transport mechanism of collisional transfer of energy. In classical fluids both mechanisms are 
present. The former is dominant at low and moderate densities and the latter at liquid densities. 

The overlapping Lorentz gas also has a percolation threshold, where diffusion only occurs below 
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differences between both modeis for percolation diffusion are discussed. The diffusion phenomena 
in both models occur on random percolation structures, which are essentially each other's duals Q. 
In Lorentz gases the return probability has a LTT ~ t' d / 2 , and the VACF has a LTT ~ t~ 1 - d / 2 . 
In the hydrodynamic interpretation of the LTT of the VACF in classical fluids in terms of vorticity 
diffusion, as given in by Alder and Wainwright, the velocity auto-correlation function of a tagged 
particle may also be interpreted as the return probability of the initial momentum to the tagged 
particle itself, where d is the spatial dimension of the system [53, 1 1 (1 Hl| . A further property of 
the Lorentz gas is that the power law tail is absent when the scatterers are placed on a periodic 
lattice [l2T | . The present paper also analyzes what happens to the long time tails when the particles 
are put on a lattice. In fact, the interpretation of C e (t), being a return probability in a diffusion 
process, already suggests that C e (t) ~ P(R = 0,t) should also have a power law tail~ t~ d / 2 on a 
regular lattice. 

As the long wave length decay modes of classical fluids and DPD models are essentially the 
same, the decay of the corresponding equilibrium correlations should also be the same. This is in 
fact implied by Onsager's regression hypothesis 01 . For the VACF in DPD fluids [l4| and for 
the energy autocorrelation function (EACF) in a DPD solid more specialized mode coupling 
arguments have been developed to show the existence of a power law tail ~ t~ d l 2 . Unfortunately 
the predictions, in particular the one for the DPD solid, could not be confirmed by the existing 
computer simulations of three-dimensional DPD systems, because the systems studied were too 
small, e.g. N = 1000, L/r c = 5 and N — 2000, L/r c = 6, where V = L d [l5|. Consequently, 
the short time kinetic decay, exp[— Wot], and the slow decay of the macroscopic diffusive modes, 
exp[—k 2 Dt], were equally fast (no separation of time scales), and power law tails in the EACF 
could not be observed because of strong finite size effects. The system sizes in our present computer 
simulations of the two-dimensional DPD solid are sufficiently large to observe the power law tails. 

The plan of the paper is as follows. In section II the definitions and results for the DPD solid 
are briefly summarized in so far as needed in the present paper. In section III the mode coupling 
theory of 0,0 is applied to obtain the explicit form of the LTT of C e (t) in the random DPD solid. 
In section IV the results are compared with computer simulations of the two-dimensional case. In 
section IV it is also explained how the same results for the time correlations can be obtained from 
deterministic simulations, where the dynamics is free of statistical noise. In section V the LTT's 
of time correlation functions are studied on a lattice. Conclusions are presented in section VI. 



II. RANDOM DPD SOLID 



The system consists of N = nV point particles at fixed random positions Yi(i = 1,2, ... , N), 
contained in a volume V = L d , and obeying periodic boundary conditions. Each DPD particle 
interacts with all particles inside its interaction sphere of radius r c . The only dynamical variables 
in the system are the internal energies ti of the particles, and the time evolution is given by the 
Langevin equation, 

3 3 

Here w(r) is a positive interaction function of finite range r c , normalized as J rfrw(r) — r d . In 
the present paper w(r) — const 9(r c — r) is proportional to the unit step function, which vanishes 
for r > r c . The kinetic rate constant Ao is a model parameter that determines the interaction 
frequency. The prime on the summation sign indicates that j ^ i. Here i*y = —Fji is Gaussian 
white noise, whose explicit form is given in paper I. It satisfies the detailed balance conditions, 
guaranteeing that the system reaches thermal equilibrium. 

To discuss the short time decay of an energy fluctuation we consider 

v = i e ( r c ~ nj)) = P^d, (2) 

3 

where v is the mean number of j-particles inside r,j < r c that interact with the i-th particle. 
Moreover Wd = it d / 2 /T{\ + d/2) is the volume of a d— dimensional unit sphere (d — 1, 2, • • •), and 
p = nr d is the reduced density. 
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The average (• • • ) denotes a quenched average over the random configurations of the fixed DPD 
particles. The basic relaxation rate at short times can be estimated from and (J2J as 



= Ao2^ (w(rij)) = pXa- (3) 



Each DPD particle is a mesoscopic subsystem with a density of internal states ~ e a , where a is 
proportional to the number of internal degrees of freedom of a DPD particle, and satisfies a> 1. 
In thermal equilibrium the single particle distribution function is 

Me) = f ^ T) m a eM-M, (4) 

where T = is the temperature. Also note that the evolution equation conserves the total 

energy, i.e. £j(t) = E = const, as both terms in £Q) are anti-symmetric in i and j. 

A final comment regarding the stochastic differential equation which contains multiplicative 
noise. As discussed in paper I, the difference between the Ito and Stratonovich interpretation of 
the Langevin equation and the corresponding Fokker Planck equation can be neglected to leading 
order in 0(1/ a) in the present model. 

In this paper we study the long time tail of the energy autocorrelation function (EACF) in 
thermal equilibrium, 

a(*) = ^£(<M*)Mo)> 0! (5) 

i 

where (• • • )o ^ s an avera g e over the canonical distribution FT . ip(tj). Here 5ej (t) = ej (t) — (ej) is the 
energy fluctuation with (ej) Q ~ a/ (3. It is essentially the return probability of the initial energy 
Sei(Q) to the particle, from which it originated. It is the analog of the velocity autocorrelation 
function, C v (t) = (v x (t)v x (0)) Q in the same sense as the VACF is the return probability of the 
initial momentum. The latter picture explains its LTT ~ i~ d / 2 [y]. However the time integral of 
the VACF equals the coefficient of self diffusion D, whereas the time integral of the EACF is not 
related to the heat diffusivity Dt, or to any other transport coefficient. 

The explicit form of the power law tail of the EACF has been given in Refs.0,^^- The derivation 
can essentially be copied from the corresponding derivation for the VACF in fluids Q. We only 
give an outline. We first express <5ei(t) = J drSe s (r, t) where 8e s (r, t) = 5ei(t)S(r — ri) is the local 
energy density of tagged particles. According to the Onsager regression hypothesis, the average 
decay of fluctuations around thermal equilibrium follows the macroscopic approach to equilibrium. 
This implies that a local fluctuation Se s rapidly decays to its value in local equilibrium, i.e. 

= j dv5e s {r,t) ~ J drn s (r)CST(r,t) ~ - ^ ST k (t)n s _ k , (6) 

k 

where C — aks is the specific heat per DPD particle, and ST(r, t) is the local temperature 
fluctuation, n s (r) = S(r — r±) the quenched local density of tagged particles, and denotes 
the Fourier transform of a(r). The subsequent slow evolution is controlled by the heat mode 
<$Tk(t) = STk(Q) exp[— k 2 DTt\. The tagged particle density n k is a static mode. Inserting 5T)~ in 
(0 and ©, and using the relations n s k — * 1 as k — > 0, and STk(0) — > 5ei(0)/C as k — * 0, we obtain 
in the long time limit, 

C e (t) = -L^ex P hfc 2 ^ T i]C e (0). (7) 

k 

This yields in the thermodynamic limit as N = nV — * oo at fixed n, the power law tail, 

C ' ' ' 1 f ^ ex P [-k*D T t] = - 1 = - L_ (f » l) (8) 



C e (0) nj {2ir) d 1 n{4TrD T t) d / 2 p(Ant*) d/2 

In the first step we used the relation C e (0) = ((<^e) 2 ) — ct/(3 2 , and in the last we introduced the 
dimensionless time t* = Dxt/r 2 , where to — f^/Drip) is the characteristic time for heat diffusion. 
The integral representation ©, / dkP k (t) — P(R = 0,t) = (47r£) T ^) _d/2 , of the long time limit of 
C e (t) shows that P(R,t) is the solution of the heat diffusion equation, and that P(R = 0,t) is a 
return probability. 

The main goal of this paper is to test the theoretical prediction about the long time tail (LTT) 
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III. LONG TIME TAILS IN STOCHASTIC AND DETERMINISTIC SIMULATIONS 

The goal of this section is to compare the predictions JJJ for the long time tails with the results 
of computer simulations in the two-dimensional random solid described by the Langevin equation 
(Q. The values of the heat diffusivity, to be used in the comparison, have been obtained from 
computer simulations, using the methods of paper I, and agree well with the mean field result 
Doo(p) — ^oP r c/[2{d + 2)] at high densities. As p decreases, -Dt(p) decreases faster than linear, 
and vanishes at a non vanishing threshold density p p , which coincides with the percolation threshold 
for continuum percolation of overlapping spheres. In two dimensions it has the value p p ~ 1.43629 
[l6j . and in three dimensions one has p p = 0.65296 [l7) . 

To measure LTT's of equilibrium time correlation functions the DPD solid has to be in a state of 
thermal equilibrium. This is achieved by initializing the system in a state with a uniform energy 
distribution ei(0) = E/N, where E is the total energy. In principle the Langevin forces, which 
satisfy the detailed balance conditions, drive the system to a state of thermal equilibrium, described 
by ipo in Eq.Q. This is only the case in a pure ergodic system. However with decreasing density 
larger fractions of particles are contained in small disconnected islands, which have no interactions 
with the bulk of the particles in the system, i.e. belong to small independent ergodic subsystems, 
and cannot redistribute their initial energies over the bulk system. Hence, or initialization of the 
system in the canonical equilibrium state can only be guaranteed if the fraction of island particles 
is negligible, i.e. for p ^> p p . 




Figure 1: Simulation data for EACF, C e (t) versus t at different reduced densities p (defined below Eq.(J5J) 
in the two-dimensional random solid for a system of N = p(L/r c ) 2 — 5 x 10 4 particles, compared with 
theoretical predictions: the dashed lines represent the short time exponential decay, and the solid lines the 
algebraic LTT~ 1/t. The higher the density the sooner the LTT is reached. 




t* = D T (p)t/r c 2 

Figure 2: Collapse plot of the LTT's in Fig.l, obtained by plotting pC e (t)/C e (0) in Eq.© versus t* , 
compared with the predicted LTT (solid line). Note that the crossover from exponential to algebraic decay 
is at to — t^/DtIp). 
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t* = D T (p)t/r c 2 

Figure 3: A more sensitive comparison of the collapse plot with the theoretical prediction (solid horizontal 
line), focusing on relatively short times. 

Figs. 1,2,3 show the simulation data for the energy auto-correlation function C e (t) of a DPD 
particle in a two-dimensional system of N = p(L/r c ) 2 particles at different densities. There is 
excellent agreement of computer simulations with theoretical predictions. In the simulations we 
measure length in units of r c , and times in units of 1/Ao. At short times the decay is exponential, 
exp[— oJot], with a rate constant uiq — pAo, as given in (jSJ), in very good agreement with the 
simulation data. At large times the plots show the long time tail ~ t~ d l 2 . 

By plotting the simulation results as pC e (t)/C e (0) versus the dimensionless time t* = DT(p)t/r 2 
the data for different p can be collapsed on a single LTT-curve. Moreover the combination of Figs.l, 
2 and 3 shows that the crossover time from exponential short time decay ~ exp[— u>ot] to power law 
decay is for all densities given by to — r 2 /DT{p). 

The finite size effects on the dynamics are controlled by the ratio L/r c = (N/p) 1 / d . As DPD 
particles in absence of conservative forces are point particles, the density can become arbitrarily 
high. So, as p increases at fixed N the finite size effects increase. Consider the data in Figs.l, 2 
and 3 at (p — 5,L = 100), corresponding to JV = 5 x 10 4 particles. If the number of particles 
is reduced to A" = 10 4 a faster decay than 1/t decay becomes noticeable, which is statistically 
significant. For N — 10 3 the decay is even faster, and looks exponential. It is caused by finite size 
effects. We return to this point later on. 

The stochastic simulations above were described by the Langevin equation Q). The subsequent 
deterministic simulations of this section correspond to the same equation with the Langevin forces 
switched off (Fij(t) = 0). Moreover, the stochastic forces have no effect on the decay of the 
equilibrium time correlation functions (5a(t)5a(0)) , provided both types of correlation functions 
are at t = in thermal equilibrium. This observation follows from Onsager's regression hypothesis 
on the average decay of fluctuations. The purpose of driving the system by Gaussian white noise, 
that satisfies the detailed balance conditions, is to maintain the system in thermal equilibrium, 
but does not affect the decay of these functions. 

The above observations offer the interesting possibility to measure the LTT of equilibrium time 
correlations deterministically, i.e. the dynamics is free of statistical noise, provided the system is 
prepared initially in the proper thermal equilibrium state. Here the thermal fluctuations are only 
accounted for in the initial distribution. This method has also been used in 0] where the kinetic 
part of the stress-stress correlation was measured using the lattice Boltzmann equation. 

Next we discuss the results of the deterministic simulations of C' e (t) in the two-dimensional ran- 
dom DPD solid at density p = 5, where on average irp particles are surrounding the central particle 
in the interaction sphere. In order to simulate the equilibrium time correlation function C e (t), we 

prepare the system in the initial state, described by the A^— particle distribution, rL=i[V'o(e*)/K|j 
and average over different runs with independent initial configurations. In spite of the additional 
average over different runs the computational effort for the deterministic simulations is consider- 
ably smaller than for the stochastic ones. Figs. 4 and 5 show log-log plots of C e (t)/C e (0) versus 
t at different system sizes L/r c with A" = p(L/r c ) d particles. Fig. 4 is focusing on the short 
time behavior exp[— Xopt] and power law tails (1/ p)[4irD T (p)t/r 2 ]- d/2 , and Fig. 5 on the power 
law tail and the ultimate exponential decay. Both figures also show for comparison the stochastic 
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and the intermediate power law behavior of the deterministic and stochastic simulations are the 
same, and agree with the theoretical prediction, until the time where the statistical uncertainty 
in the stochastic simulations has become too large. This behavior beautifully confirms Onsager's 
regression hypothesis. 




Figure 4: The log-log plot shows the EACF vs t in the deterministic simulations at p = 5 for L/r c — 
20, 40, 100, corresponding to N = 2 x 10 3 , 8 x 10 3 , 5 x 10 4 particles. The scattered black circles (L/r c = 44) 
show stochastic simulations for comparison. The straight solid line shows the predicted power law tail JSj, 
and the dashed line indicates the short time exponential decay. 




Figure 5: The log-log plot shows the same deterministic simulations as in Fig. 4, but focuses on the 
intermediate algebraic tail. The straight line is the LTT in the thermodynamic limit. The ultimate long 
time decay is again exponential. For the scattered black circles see caption Fig.4. 



Fig. 6 shows that the ultimate decay is again exponential. This is a finite size effect, which 
is in excellent agreement with the mode coupling prediction (JJJ for finite systems, where the 
k-summation can not be replaced by a k-integral. On a square lattice with periodic boundary 
conditions k a — 2im a /L with a — x,y and n a = 0, ±1,±2, The integral approximation 
to (J7J is only valid in the limit of large t, and small k, such that k 2 t — constant. This asymp- 
totic approximation is appropriate as long as Dxk^^t -C 1, and breaks down at t m in f° r ; say, 
Dxk^^tmin m 1/3, where the minimum fc-value is k m i n — 2ir/L. The first few terms in the finite 
size mode coupling formula on a square lattice follow from J7J in the form, 

tf| = jf E -Ph |n| 2 S ] = 1 (4e- + 4e- 2 * + 4e^ + 8e^ + ■■■), (9) 

where s = k^^Drt — 4TT 2 Drt/L 2 , and n = (n x ,n y ). In Fig. 6 the straight lines of simulation 
data are well represented through the first term (4/iV)e~ s of the right hand side of @. 

Furthermore, our method is not suitable to study equilibrium time correlations close to the per- 
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Figure 6: The half- logarithmic plot focuses on the ultimate exponential decay in the simulations of Figs. 
4 and 5, which are finite size effects, and are quantitatively accounted for by the first term in Eq.@, 
(4/A) exp[— (2n/ L) 2 Drt]. The solid line at the top right is again the algebraic LTT for thermodynamically 
large systems. For the scattered black circles, see caption Fig.4. 

the finite fraction of particles contained in smaller disconnected islands form uncoupled ergodic 
subsystems. This prevents the system from reaching the canonical equilibrium state, • tpo(ej). 

IV. DPD - SOLID ON A LATTICE 

In this section we introduce a lattice version of the DPD solid in QJ, where N particles fill the 
N sites of a hypercubic lattice with lattice distance a, and volume V = L d = Na d . In this model 
the interaction range r c is the control parameter and heat diffusion only occurs for r c > a. We 
will observe that the LTT of C e (t) in the random DPD solid survives when the DPD particles are 
put on a periodic lattice, as is expected for a return probability. To do so we consider the special 
case of the DPD solid on a lattice, where C e (t) can be calculated exactly la. H owever, in periodic 
Lorentz gases the power law tail disappears in the VACF, as shown in Ref.|l2j|. 

Consider the evolution equation for the lattice space-time correlation function C(n, t) = 
(<5 e„(t)<5 cq (0))o in thermal equilibrium, where n = (n x ,n y , . . . , no) labels the sites of a hypercubic 
lattice with periodic boundary conditions, and n a — {0, 1, 2, . . . , L — 1} with a = {x, y, . . . , d}. 
The evolution equations are obtained by replacing e; in (|TJ with Se n (t), multiplying that equation 
with <5eo(0), and averaging over the iV-particle equilibrium distribution function. As <5en(0) is 
uncorrelated with the Langevin force, the equations of motion of the lattice space-time correlation 
functions follow as, 

dC(n, t)/dt = A ^(r nm ) (C(m, t) - C(n, t)) (10) 

m 

with initial condition C(n, 0) = <5 n o((<5e) 2 )o ~ 8 n oa/(3 2 . This is a discrete diffusion equation, 
which can be solved exactly for general dimensionality [5j , and the single site correlation or return 
probability C(n = 0,t) = C e {t). For w(r nm ) restricted to nearest neighbor interactions (r c = a), 
the solution is given by Q , 

C e (i)/C e (0) - [e- T I (r)] d - (2^r)- d / 2 . (II) 

Here — * denotes the long time behavior for t> 1. Moreover Ii(t) with I = is a modified Bessel 
function with r = 2DTt/a 2 , and the heat diffusivity, Dt = Xo<X jvjfa where units of time are 
consistent with and ©. 

In the lattice model both characteristic time scales, the short kinetic time scale, 1/ujq = l/pAo, 
and the long diffusive time scale to — r c/^T(p), are roughly inversely proportional to the mean 
number of j-particles inside the interaction sphere (r.y < r c ) of the central particle i, and dif- 
fer roughly by an order of magnitude in size As p J, p p this number decreases rapidly, and the 
cauilibration time tn at o„ diverges, as the interactions become weaker and weaker. 
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For the exactly solved case of the lattice DPD solid in d-dimensions with nearest neighbor 
interactions, Figs. 7 and 8 show for the one-dimensional case that a very long equilibration time 
Aoio — 5 x 10 5 is required. After that period the exact solution C e (t) can be observed with its 
full tail over the whole time interval Xot < 1000. If the equilibration time is too short, the 

observed LTT falls off faster than 1/yt . 




Figure 7: The plot shows the simulated EACF on a one-dimensional lattice with nearest neighbor interac- 
tions after different equilibration times to — 5, 50, 500, 10 4 and 5 x 10 s , as well as the exact solution 11 U 
for the one-dimensional lattice model. 
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Figure 8: The plot shows the same data for C e (t) as in Fig. 7, but focuses on short times. Note that the 
interval around Xot ~ 1 also equilibrates very slowly. 



For general interaction ranges (r c = Ma), where M represents the number of interacting lat- 
tice sites on periodic lattices, the mode couplings arguments of @, produce the LTT in JSJl 
which yields in one dimension an algebraic LTT ~ \j\fi. In the one-dimensional lattice case it is 
straightforward to calculate the coefficient of heat diffusivity Dr{p) analytically, as shown in the 
appendix, and the results of the LTT's in the stochastic simulations are show in Figs. 9, 10, 11, 
which can be extended directly to higher dimensional lattices. The simulation results are in very 
good agreement with the predictions of the mode coupling theory. 

We also note that the random DPD model in one dimension is in the thermodynamic limit 
nonconducting because at any given density p there will always be a iron- vanishing probability to 
find a single pair (ij) of nearest neighbors with > r c . 



V. CONCLUSIONS AND OPEN PROBLEMS 



The general conclusion from the previous sections is that the predictions of mode coupling theory 
for classical fluids regarding the existence of LTT ~ t~ d / 2 in time correlations, are in excellent 
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Figure 9: The EACF C e (t) in the one-dimensional lattice model at different p for a system of N = 10 4 
particles compared with the theoretical predictions: dashed line is the short time exponential decay, and 
the straight lines are the algebraic LTT ~ 1/vt (Compare with Fig.l). 
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t* = D T (p)t/r c 2 

Figure 10: Collapse plot of the LTT in the one-dimensional lattice model (Compare with Fig. 2). 



function C e (t) = (8ei(i)5ei(0))o in the random DPD solid and in the lattice DPD models (with a 
uniform density distribution) both in one and two dimensions. Here the energy auto-correlation 
function C e (t) at large time is essentially a return probability in the same sense as the velocity 
auto-correlation function is in a classical fluid. 

Nevertheless, a number of interesting questions about consistency of the theory remains for 
the LTT-predictions C e (t) oc (Z>rt)~ d / 2 for d — 1,2, ... . For classical fluids in two dimensions the 
Navier Stokes transport coefficients do not exist. For example, let C(t) — (v x (t)v x (0))o with a LTT 
~ f -d / 2 , be the velocity autocorrelation function. Then the long time limit of D(t) = f Q cLtC{t) 
yields the coefficient of self diffusion D, if the integral exists. However, this relation would lead 
to the result D(i) ~ {V~t, \nt} for d = {1,2} respectively in the long time limit, leading to the 
non-existence of the Navier- Stokes transport coefficients in one-dimensional and two-dimensional 
fluids. Then, self consistent mode coupling theory, ring kinetic theory and other renormalization 
procedures l e& d to a "renormalized super-long time tail" of the form pH 12(11] . 

D(t) ~ {t 1 ^ 3 ,^} 

C{t) ~ dD(t)/dt = {t- 2/3 ,l/[tVhTi}} (12) 

for d = {1,2} respectively. In fact, van der Hoef and Frenkel [23] have confirmed for the VACF 
in two-dimensional fluids the existence of a "faster-than-i _1 -tail" of C(t) by computer simulations 
using lattice gas cellular automata. However, our computer simulations of C e (t) in Figs. 1-3 and 
9-11 do not show any deviations from the predictions of the simple mode coupling theory in 
and seem to confirm the LTT ~ t~ d / 2 , as well as the finiteness of Dt- 

This simple behavior is further confirmed by the observations (i) that C e {t) in the two- 
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Figure 11: A more sensitive comparison of the collapse plot of Fig. 10 for short times (compare with Fig. 3) 
with the theoretical predictions. 



10, 11 behave essentially the same, and (ii) that the lattice models contain as a special case the 
model with nearest neighbor interactions (with r c = a), which is exactly soluble for all values of 
d 5], and show the same LTT oc [DTt)~ d l 2 with a finite Dt, as the random DPD solid and the 
remaining lattice models with r c = Ma, where M — 1,2,.... In fact, the behavior of the DPD 
solid resembles that of Lorentz gases, where the diffusivity D is also finite for all values of d. In 
fact the random DPD solid and the Lorentz gas are in several respects each others duals , where 
the VACF has the LTT - t" 1 "'*/ 2 . 

Furthermore, there is no contradiction between the long time behavior of fluids in (|12l) and 
that of the DPD solid in (JSJ), because Dt in the latter case is not related to the time integral of 
C e (t), but is given by the Green-Kubo formula, involving the time correlation function Cq(£) of 
the microscopic heat current, whose time integral converges |f§- 
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Appendix A 

In this appendix we calculate the coefficient of the heat diffusion for the lattice version of the DPD 
model in one dimension with an interaction range r c = Ma with M = 1, 2, • • • for the evolution 
equation with the Langevin force set equal zero, and using the range function w(r) ~ \9{r c — r) 
with [w] = r c for the one-dimensional case. The evolution equation JQ) for r nm > Ma becomes, 

3e(r n . 1 ^ — r . ^ 

— Q- t — =2 ° [e[r m +ma,t) - e(r n ,t)\, (A.l) 

\m\<M 

where M — ±1 for nearest neighbors. For large spatial scales, this discrete version of the diffusion 
equation is expanded to C(V 2 )-terms included, with the result 

^WvMM) (A.2) 

with a heat diffusivity 

1 M 1 

D T = ttAqo 2 V m 2 = —\ a 2 M(M + 1)(2M + 1) (A.3) 
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It is convenient to eliminate a and M from this expression in favor of the interaction range r c = Ma 
and of the reduced density p. Here pzu^ (wi — 2) is according to © the number of particles 
interacting with a given particle, which is here 2M. So, p = M, and the heat diffusivity becomes, 

^)4pW(l + I)(l + l). (A.4) 

In the limit of large density the heat diffusivity approaches the value D oc (p) = ^p\ r^. In Fig. 12 
the analytic result for R\(p) = Dt(p)/D oc (p) = (1 + l/p)(l + l/2p) is compared with computer 
simulations of the DPD-solid in a one dimensional lattice and shows excellent agreement. We note 
that for the case of n.n. interactions (p = 1) the heat diffusion becomes Dr(p = 1) = ^XqO 2 . 



3 - 




Q I i i i i I i i i i I i i i i I i i i i I 

5 10 15 20 

P 

Figure 12: A comparison of the simulations with the theoretical predictions for the heat diffusivity in the 
one- dimensional lattice model for Ri(p) = Dt (p) /-Doo (p) = (l+-)(l + j^) with a few simulation data. 
Note that p is an integer. There is excellent agreement. 
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